Dipolar Dynamics for Interacting Ultracold Fermions in a Trapped Optical Lattice 
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Using the time-dependent density matrix renormalization group (tDMRG) method, we calculate transport 
properties of an interacting Fermi gas in an optical lattice with a confining trap after a sudden displacement 
of the trap center. In the regime of attractive interactions, the dipolar motion after the displacement can be 
classified into underdamped oscillations and overdamped relaxations, depending on the interaction strength. 
These numerical calculations are consistent with experimental results. In the regime of repulsive interactions, 
we predict a revival of the oscillations of the center of mass when the interaction strength is increased. This 
unique feature can be considered as a dynamical signature for the emergence of a Mott plateau for an interacting 
trapped Fermi gas in an optical lattice. 
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I. INTRODUCTION 

Ultracold atoms have proven to be an ideal platform to 
study a number of unsolved problems in quantum many-body 
physics. The unprecedented controllability of the system is 
experimentally achieved through optical lattices [1] and Fes- 
hbach resonances [2]. Optical lattices provide the possibility 
to manipulate the dimensionality of the lattice and the ratio of 
interaction to kinetic energy, while Feshbach resonances can 
be used to tune the interaction strength, crucial to many-body 
effects. All these developments provide a route to simulat- 
ing many paradigmatic quantum models in strongly correlated 
systems and therefore are a novel approach for answering fun- 
damental problems in many -body physics [3]. 

Driven by the ambitious goal of carrying out quantum sim- 
ulations, experimentalists have begun to explore the underly- 
ing physics of strongly correlated fermions in the atomic ap- 
proach. The Hubbard model is a simple but very good start- 
ing point for modeling the essence of interacting fermionic 
systems. In the meantime, great success has been achieved to- 
wards experimental realization of this model. The Fermi sur- 
face of degenerate Fermi gases has been observed in a three- 
dimensional optical lattice [4]. Molecules from fermionic 
atoms have been created with a Feshbach resonance [5], and a 
Mott insulator has been detected [6, 7]. Although all these de- 
velopments have paved the way for understanding fermionic 
superfluidity, transport properties are hard to access. 

In one approach to transport phenomena, Strohmaier et 
al. [8] carried out an experiment on an interacting Fermi gas 
in an optical lattice. They used the center of mass (COM) mo- 
tion of an interacting Fermi gas to characterize its dynamical 
features. More specifically, the COM motion of the gases was 
monitored after a sudden shift of the external harmonic trap 
by a few lattice sites, giving rise to dipole excitations. Unlike 
similar experimental studies on bosonic atoms with repulsive 
interactions [9], they focused on the regime of attractive inter- 
actions and found that with increasingly attractive interactions 
weakly damped oscillation turns into a slow relaxational drift. 



A quantitative understanding of all these observations is a 
theoretical challenge due to the non-equilibrium nature of the 
process. In bosonic experiments [9], the dipolar motion due 
to the sudden shift reflects the dynamical excitations in the 
many-body system, and has stimulated considerable theoret- 
ical interest, involving semiclassical solutions [10] and nu- 
merical diagonalization of the Bose-Hubbard Hamiltonian in 
a small system [11]. By using the time-dependent density- 
matrix renormalization group technique (tDMRG) [12-17], 
simulations of the bosonic model have been successful in 
quantitative comparisons with the experimental data [18, 19]. 
This indicates that tDMRG is an ideal tool in the strongly cor- 
related regime in one spatial dimension even far away from 
equilibrium. 

In this work, we present a comprehensive numerical sim- 
ulations of the fermionic experiment by Strohmaier et al. [8] 
in the attractive interaction regime, and extend the study to 
the repulsive interaction regime, which is also experimentally 
achievable. To take into account the full time dependence, we 
focus on one spatial dimension by using exact diagonalization 
and tDMRG. Our approach is similar to that of Okumura et 
al. [20] who present results for the same system, but we go far 
beyond this short note by analyzing the results and explaining 
the qualitative and quantitative features and studying the re- 
pulsive case. In our tDMRG simulations, up to 400 states are 
kept in the reduced Hilbert space, and a second order Trotter 
decomposition is used in the time evolution. 



II. MODEL 

Ultracold fermions in a one-dimensional optical lattice can 
be described by the Hubbard model [21, 22] 

H{t) = -J I>l+i,<Aa + Hc ) + u H "<.t"a 

+ VL^h i ^[i-i (i + Q(t)&x] 2 . (1) 



2 



The first term describes the tunneling of fermions between 
nearest neighboring sites, where J denotes the tunneling 
matrix element between adjacent lattice sites and Cj iCT the 
fermionic annihilation operator on site i in the "spin" state 
a (f or I). The second term is the on-site interaction with 
strength U, where ni jCr = c| cr c j a . The last term models the 
additional confinement of the harmonic trap with curvature f2, 
initial center in and displacement of the trap center Ax. In all 
our calculations, the shift Ax is set to be 5, close to that in the 
experiment [8]. For brevity we work in units where the lattice 
spacing is unity hereafter. To study the dipolar motion, we are 
interested in the position of the COM, which is defined as 
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Here Nf = iVj, = TV is the particle number per "spin" species 
of a population balanced gas. 

III. ATTRACTIVE INTERACTIONS 
A. Weakly interacting regime: underdamped oscillation 

Unlike the bosonic case, the fermionic dipole oscillations 
are weakly damped even in the non-interacting case. This 
damping is caused by the dephasing of particles in the same 
"spin" state due to Pauli's principle. When a weak on-site 
interaction is introduced, the dephasing will be more pro- 
nounced because interaction between particles with different 
spins also contributes to incoherence. 

Fig. 1 (a) shows these oscillations in a system with L = 80 
and N= 16. In this regime the COM motion can be described 
by an underdamped harmonic oscillation around the new trap 
center after the shift. We fit the motion of COM to 
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FIG. 1. (Color online) (a) The position of COM as a function of 
time with different interactions, (b) U / J dependence of the damping 
rate /3. The line is a fit of the data to Eq. 4. The trap curvature is 
fi = 0.005J, equivalent to a trapping frequency of about 50 Hz in 
the experiment. 



X(t) = A + Be~ 0t cosu;t 



(3) 



where /? is the damping rate and lu is the frequency. Fig. 1 
(b) shows the damping rate f3 as a function of the interac- 
tion strength U/J. The results are well fitted by the empirical 
power law 
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with a = 0.0104, b = 1.62, and c = 0.00258. The first term is 
a power law describing the interaction effect and the second 
term describes the dephasing due to Pauli's principle. 



B. Strongly interacting regime: overdamped relaxation 



this critical value is estimated by examining whether the COM 
oscillates around the new trap center or not after the shift. 

As can be seen in Fig. 2 (a), the motion of the COM in 
the strong interacting regime is characterized by a relaxation 
towards the new trap center. Unlike the previous weakly inter- 
acting regime, the COM mass can no longer oscillate around 
the new trap center, although its motion still shows some neg- 
ligible undulations (see the case U — —8J in Fig. 2 (a) for ex- 
ample). As the attraction increases the fermions form stronger 
bound pairs with a larger effective mass. This slows the mo- 
tion of the atoms towards the new trap center. 

We extract the relaxation rate by fitting the curves to an 
exponential decay 



z(t) = Zoo (l-e- rt ), 



(5) 



In the strongly attractive case the experiment [8] shows a 
relaxational COM motion. Here we study this overdamped 
behavior in a system with L = 80 and il = 0.005 J. First, we 
find a critical value U c « —3.5 J for the crossover between 
underdamped oscillations and overdamped relaxation in sim- 
ulations with 7V=10-20 particles per "spin" species. Note that 



where z,^ is the position at t oo, which is assumed to be 
the new trap center (in our case, = —5). 

The relaxation rate as a function of interaction for different 
particle numbers per "spin" is shown in Fig. 2 (b). We find 
very good agreement with experimental data in the regime of 
strong interactions. 7V=16 corresponds to a system with half 
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FIG. 2. (Color online) (a) COM position as a function of time in a 
system with N = 16. The interaction strengths U are indicated in the 
figure, (b) The relaxation rate F vs the ratio U/J for various particle 
numbers per "spin" in the lattice. The experimental (black circle) 
results are obtained from Fig. 3 of Ref. [8]. Note that our results can 
be qualitatively compared with the experimental data although the 
experiment was performed in a three-dimensional optical lattice. 



filling at the center, consistent with the experiment. We also 
find that for a given interaction strength U, the relaxation rate 
decreases as the particle number increases, as has been ob- 
served experimentally. The good agreement is a strong evi- 
dence that the single-band Hubbard model can capture all the 
features of the experiment. 

In the strong interaction limit the fermions form local sin- 
glet pairs on a lattice site. The transport properties in this 
regime are thus governed by the dynamics of these local pairs, 
which can be approximated as hard-core bosons (HCB). Here 
we perform simulations to test the validity of this mapping to 
an effective HCB model 

+ O'^Sj6 ?; ( l -z -e(i)Ax) 2 , (6) 

i 

where &i's are bosonic annihilation operators with additional 
hard-core constrains b\ 2 =Sf = and {b i , b\} — 1. The effec- 
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which reduces to 2J 2 /\U\ in the large \U\ limit [23]. Since 
one HCB represents a pair of fermions, whose bare mass is 
twice that of a fermion, the trapping curvature ft' = 2S7 is 
twice that of the original one. 

By exactly diagonalizing the model it is straightforward to 
obtain numerically exact results for the time evolution of the 
system. 
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FIG. 3. (Color online) Static and dynamical mappings to the HCB 
model. Density distributions in the initial states for (a) U = —24 J 
and (b) U — — 5 J, respectively, (c) The relaxation rate T as a func- 
tion of U/J. The results of both the original fermionic model and 
the HCB model are shown for comparison. Here we use N — 16 
particles per "spin" species. 

To test the validity of the HCB approximation for static 
properties, we consider the density distribution. From the 
comparison of the fermionic and the effective model shown 
in Fig. 3 (a) and (b), we find that the mapping of static quan- 
tities to an HCB model is valid in the large \U\ limit, but less 
so in the regime of intermediate interactions. The HCB model 
overestimates the central density at U ^ — 5 J because it as- 
sumes most fermions form fully local pairs even when the at- 
tractive interaction is not strong enough. For the dynamical 
mapping, the HCB model captures the interaction dependence 
of r qualitatively. However, the poor mapping of static quan- 
tities at U ~— 5 J significantly compromises the approximate 
relaxation rate T based on the HCB model, which deviates 
from the exact results by about 50%. The HCB model also 
underestimates the relaxation rate in the intermediate regime, 
where single particles and local pairs coexist. 
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We can improve the approximation by using an extended 
HCB Hamiltonian 

H' B (t)=H B (t) + V^n i f k+1 , (8) 

i 

where the term with V > describes the nearest neighbor in- 
teraction between local pairs. This repulsive term will reduce 
the deviation in density profiles shown in Fig. 3 (b). How- 
ever, the dynamical mapping cannot be perfect even with this 
extended HCB model, because the internal dynamics of the 
pairs are not taken into account. 

IV. REPULSIVE INTERACTIONS 

In the repulsive case we need to distinguish between low 
and high density regimes. At low density (or in a shallow trap 
at fixed particle number), the fermions are delocalized, lead- 
ing to compressible metallic states with central density less 
than unity at U = 0. In this case, Mott-insulating states will 
never appear even when U is increased due to the low density. 
On the other hand, if the particle density is high enough (or 
the external trapping potential deep enough), a central Mott 
domain will be formed beyond some critical value of U. 

In the low density (shallow trap) case, the COM motion is 
always underdamped, as shown in Fig. 4 (a) and converges in 
the large U limit. 

However, for a high density (deep trap), the situation is very 
different, as shown in Fig. 4 (b). Here, there is a crossover 
between underdamped and overdamped behaviour at around 
U = 2 J. However, when U is further increased beyond 4 J, a 
Mott domain will be formed at the center (see Fig. 5 for ex- 
ample), and oscillations reappear on the background of a very 
slow relaxation. By increasing U the dipolar motion exhibits 
a collapse and revival of oscillation with interaction strength. 

Analyzing the revival of the COM oscillations for SI = 
0.01J and U = 12J as shown in Fig. 5, we see that their 
frequency is the same as the non-interacting one, i.e., the har- 
monic trap, similar to the case of cyclotron motion of elec- 
trons [24]. Furthermore, the oscillating amplitude is about 
one quarter of the non-interacting case, suggesting that only 
one quarter of the atoms are oscillating. This agrees with the 
corresponding density profile (see Fig. 5), which shows that 
half of the atoms form the Mott plateau, leading to one quarter 
at the each edge of the domain. We argue that the whole COM 
motion can be interpreted within a two-fluid model as follows: 
The metallic component at the edges of the Mott domain oscil- 
lates with moderate damping, while the incompressible Mott- 
insulating counterpart becomes nearly localized. In this two- 
fluid model, the oscillation due to the metallic component is 
also conceptually equivalent to the residual current in the ef- 
fective model with two coupled bands separated by U [19]. 

V. CONCLUSION 

Using tDMRG, we have presented an accurate analysis of 
the dipolar motion of an interacting Fermi gas in an optical 
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FIG. 4. (Color online) COM motions for (a) SI = 0.005J and (b) 
SI — 0.01J, respectively. In each panel, the on-site interactions are 
indicated in the figure. Other parameters are: L — 80 and TV = 16. 



lattice. In the regime of attractive interactions, the numerical 
calculations for the COM dynamics after a sudden displace- 
ment of the trap minimum are in good agreement with exper- 
imental results. 

To further theoretical understanding, we have mapped the 
fermionic model to an effective HCB model in the strong in- 
teraction limit. To improve the validity of the mapping for 
intermediate interactions, an extended HCB model is neces- 
sary and will be subject of further studies. 

For repulsive interactions a very different behaviour has 
been found. In a deep trap, we find a revival of the oscilla- 
tion for the COM by increasing the interaction strength when 
the trap is deep enough to support a central Mott domain with 
increasing U. This revival is due to the fact that as the cen- 
tral Mott plateau is nearly frozen, there will inevitably mo- 
bile metallic states at the edges. By comparing with the lower 
density case, we conclude that the overdamped relaxation in 
the repulsive Hubbard model can be regarded as a dynamical 
signature for the emergence of the Mott domain for an inho- 
mogeneous Fermi gas in an optical lattice. 
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FIG. 5. (Color online) Density distribution for U — 12J. The par- 
ticle number in the Mott domain is about 16, and there are about 8 
fermions in each metallic domain on the left and right. Therefore, 
due to the imcompressibility of the central Mott phase, only a quar- 
ter of all the particles contribute to the oscillation after the sudden 
shift. 
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